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ABSTRACT 

Spectral methods have been successfully applied to the simulation of slow 
transients in gas transportation networks. Implicit time advancing techniques 
are naturally suggested by the nature of the problem. 

The aim of this paper is to clarify the correct treatment of the boundary 
conditions in order to avoid any stability restriction originated by the 
boundaries. The Beam and Warming and the Lerat schemes are unconditionally 
linearly stable when used with a Chebyshev pseudospectral method. Engineering 
accuracy for a gas transportation problem is achieved at Courant numbers up to 
100 . 
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1 * INTRODUCTION 


Spectral methods have been recently applied to the numerical simulation 
of the unsteady flow of a gas in long distance transportation networks (see 
[1]). The regularity properties of the solution [10] make spectral methods 
particularly effective for this class of problems* Engineering accuracy is 
achieved at a lower computational cost than by using more conventional finite- 
order methods. 

Transients occurring in the normal operation of pipeline networks are 
usually slow. This happens because variations imposed on the physical 
variables at the boundaries are very slow. Moreover, their propagation toward 
the interior is damped by the presence of a strong friction effect. 

These reasons, as well as the need to avoid the severe stability 
conditions which arise from the use of explicit methods for spectral methods, 
make it highly desirable to use implicit methods in time. If a method is 
considered which is unconditionally stable for the pure initial value problem, 
the boundary conditions have to be incorporated into the numerical scheme in 
such a way that no spurious stability restriction is introduced. 

Theoretical and experimental results on the numerical treatment of 
boundary conditions for finite difference approximations of hyperbolic systems 
are widely available in the literature. In [6], Gottlieb, Gunzburger, and 
Turkel address the issue of correctly imposing the boundary conditions in 
terms of the physical variables, rather than in terras of the Characteristic" 
ones. They show that the boundary conditions can be properly imposed within a 
finite difference or finite element method which is explicit in time, by a 
procedure of boundary corrections at the end of each time step. There is 


computational evidence that their precedure works for spectral methods as well 
(see, e.g., [1]). 
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It is the present authors' impression that implicit time advancing 
methods have not received sufficient attention in [6], hence the reading of 
that paper may lead to some misunderstanding about the practical 
implementation of the boundary conditions in implicit methods. As a matter of 
fact, the boundary corrections proposed there destroy unconditional stability, 
as it will be documented below. 

The purpose of this paper is to discuss in more detail the correct 
implementation of the boundary conditions within an implicit time advancing 
scheme when a Chebyshev collocation method is used in space. We stress that 
the only unconditionally stable treatment of the boundary conditions consists 
of imposing at each endpoint the prescribed physical conditions together with 
certain linear combinations of the physical differential equations. The 
coefficients of these combinations are those which express the incoming 
characteristic variables in terras of the physical ones. Thus the equations at 
the boundaries have to be incorporated into the matrix to be inverted at each 
time step. 

The Beam and Warming scheme and a class of Lerat-type schemes have been 
chosen in the following discussion as time-marching methods for the time 
discretization of a linear 2x2 hyperbolic system. An application to a gas 
transient simulation of industrial interest is also presented. 


2« THE BOUNDARY TREATMENT ON A LINEAR PROBLEM 

The following simple hyperbolic system has been widely considered in the 


literature as a model for more complex situations: 
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w t + Aw^ =0 -1 <x< 1, t>0, 


W = (Wj,w 2 ) ; A = 


l /2 1 
1 1/2 


( 2 . 2 ) 


The system is supplemented by an initial condition 


w(x,0) = w (x), -I < x < 1, 


(2.3) 


and the boundary conditions 


WjC-l.t) = WjCl.t) =0, t > 0. 


The initial-boundary value problem (2.1) - (2.4) is well-posed (see, e.g.. 


[4]); more precisely, for each t > 0 one has 


/ |w(x,t)| 2 d x <^ / |w^(x)| 2 dx. 
-1 -1 


The matrix A has two eigenvalues of opposite sign 


x = 3 x = - 1 

A 1 7 » a 2 7 * 


It can be diagonalized as 


A = TAT 


A = diag(X , ,X 9 ); T = T _1 = — ( 

1 2 / 2\1 
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The system (2.1) can be written in diagonal form as 


z t + A z x = 0 (2.9) 

where 

z = (u,v) T = T 1 w (2.10) 

are the characteristic variables. The boundary conditions (2.4) become 

(u + v)(-l,t) = (u + v)(+l,t) =0, t > 0. (2.11) 


We consider a spatial approximation of (2.1) based on the Chebyshev 
collocation method at the points 


J7T 

X. = COS ^-rr , 

J N ’ 


j = 0, • • • ,N. 


( 2 . 12 ) 


For each t > 0, the approximate solution, which we still denote by w(x,t), 
will be a couple of polynomials of degree N in the x variable. The 
degree N will be kept fixed throughout the paper. 

We denote by 


D f d ij}o<i,j<N 


(2.13) 


the matrix of the Chebyshev pseudospectral derivative at the points (2.12) 
(see, e.g. , [7]). We recall that if an N— degree polynomial w is identified 
with the vector w of its values at the points (2.12), then Dw is the 
vector of the values of w x at the same points. 
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We will now introduce time discretization methods for the previous 
system. Hereafter, we will say that a method is stable for a given At if 
the computed solutions w n at the times t = nAt satisfy an estimate of the 
form 

l!w n il 9 _< C n=l,2,3,*«* 

L Z (-1,1) 


with C independent of n. Unconditional stability will mean stability for 
all At > 0. 


We first give numerical evidence to the fact that if the boundary 
correction procedure described in [6] is applied within an implicit time 
advancing method, a severe stability restriction occurs. 

The differential system (2.1) is collocated at all the points (2.12) 
(including the boundary points) and advanced in time by one step of the Beam 
and Warming method [2], which here reduces to the classical Cranck-Nicolson 
scheme : 


~n+l 

w. 

J 


.At A 

+ -r A K 


n 

J . = w. 
J 



j = (),••• ,N. 


(2.14) 


(Let us recall for further reference that the Beam and Warming scheme applied 
to the conservation law w t + f(w) x = 0 is, before space discretization: 

Aw n + (A(w n )Aw n ) x = -Atf(w n ) x (2.13) 

where Aw 11 = w n+ ^ - w 11 and A(w) = f^(w)). 

The interior values are retained, i.e., = f° r J = 1>* ## >N-1, 

while the boundary values Wq + ^ and w^ + ^ are corrected by solving the 

linear systems 
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(wi) 0 n+1 " 0 


. n+1 


, n+1 


<*1^ + (w 2>r - ^iC + <* 2 >S 


s n+l 


( 2 . 1 ) 


and 


<“i >r - o 


'i n+1 - \ n+1 - /-~\n+l xn+1 

(W 1 } N (W 2 } N ~ (W) N " (w 2 } N 


Numerical experiments show that this method is stable only if 


At < .74 At. 


( 2 . 1 ) 


1 6 

where At* = — j is the stability limit for the Modified Euler method. 

3N 

Although the Cranck-Nicolson method is A-stable for a scalar equation, the 
stability condition introduced by the boundary treatment is more severe than 
that of an explicit method. 

The proper way to treat the boundary conditions is naturally derived 
using the characteristic form (2.9) of the hyperbolic system. In this case it 

is possible to represent the spatially discretized system as an O.D.E. system 
in the form 


dz_ 

dt + “5. " °> 


(2.18) 


where the matrix M already takes into account the boundary conditions. It 
is clear that if M is diagonalizable and has the spectrum in the right 
complex half plane, then any A-stable time discretization method yields an 
unconditionally stable scheme for (2.18). In turn, this gives rise to an 
unconditionally stable scheme for the physical system (2.1) - (2.4), via the 
transformation (2.10). 
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It is well-known that the Chebyshev collocation method at the points 
(2.12) for the scalar initial boundary value problem 

u+u=0 -1 <x< 1, t>0 

t x 

u(-l,t) = 0 (2.19) 

u(x,0) = Ug(x) 

gives stable numerical results. (Unfortunately, so far no proof of an 
estimate of the type (2.5) is known for this method. A stability result 
involving global x - t norms can be found in [8].) The method can be 
written as an O.D.E. system in the form 

dii 

— + Du = 0 (2.20) 

dt — 

T 

where u_ = ( u q >• • • ,u N _^) and ® is t * ie submatrix of D obtained by 
deleting the last row and column. Thus the differential equation is 
collocated at the interior points and at the outflow boundary point x = 1. 
The matrix D has distinct eigenvalues with strictly positive real parts, 
hence, any A-stable time discretization method will produce stable solutions 
with no stability restriction. 

In analogy with the scalar case, the spatial discretization of the system 


(2.9) is as follows: 
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The eigenvlaues of M have been computed for increasing values of N. They 
were found to be distinct and with non-negative real parts. 
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It follows that the Beam and Warming scheme will produce unconditionally 
stable approximate solutions of (2.21). They are defined by the system 


n+1 3 At f n+K n 3 At r ns 

u j + 2-(" x h -“j - nK)j 


j = 0,...,N-1 


1 it f n+ls 
)j 


n , 1 At / n% 

Vn 


n+1 _ n+1 

V 0 u 0 


j = 1,...,N 


(2.24) 


We now go back to the physical equations using the transformation (2.10). At 
each interior point both the characteristic equations are collocated, hence we 


n+1 . At .r n+K n At . / n-> 

w. + — 7T A( w 1 . = w. - — T Afw I . 

J 2 x J 2 j 2 ^ x-'j 


j = 1 , • • • ,N-1 . (2.25) 


At the boundary point x Q , equation (2.24.1) with j = 0 becomes 


t, , w + 


n+1 .At . n+1 


n+1 . At . n+1 


» i , r irrj - , at. , 
Aw x Jl + T 12 t W + ~2 Aw : 


_ r n At . tl«i . r lA nL 4 iiq 
T llL w “ T Aw U + t 12 L w ““T Aw x^2» 


x J2 


n At . n n 


(2.26) 


where we have set T 


[ T hkh<h,k<2* Similarl y at * N Eq. (2.24.2) with 


j = N gives 
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(2.27) 


(2.28) 


We conclude that at each boundary point the equation to be added to the 
prescribed boundary condition is obtained by collocating a linear combination 
of the physical equations after their discretization by the Beam and Warming 
scheme. The coefficients of this linear combination are entries of the 
matrix T”* defined by (2.10). 


A class of time discretization schemes of implicit type has been recently 
proposed by Lerat for finite difference approximations of hyperbolic systems 
(see [9], [3]). The interest of such methods lies in the fact that 
unconditional stability (in time) is achieved by including a second derivative 
term (in space) into the scheme. Unlike the Beam and Warming method, Lerat^s 
schemes are dissipative in the sense of Kreiss, hence spurious oscillations 
are automatically damped. 

We will introduce hereafter the analogs of a subclass of Lerat schemes in 
the case of spatial discretizations by the Chebyshev collocation method. The 
previous considerations on the treatment of the boundary conditions will guide 
us in defining a scheme which is unconditionally stable in time. 
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Let us recall that for the conservation law 


w t + f(w) = 0 


(2.29) 


a class of methods of Lerat's type reads as follows ([9]) 


Aw 11 + g (A^(w n )»Aw”) = -Atf(w n ) + -^1— (A^(w n )w n ) (2.30) 

“ XX X Z XX 


where Aw 11 = w n+1 = w 11 , A(w) = f'(w) and g is a negative parameter. The 
method is obtained from the two-term Taylor formula for w at t = t n . Time 
derivatives are replaced by space derivatives using (2.29) and w n is 
replaced by w n + gAw n in the second-order term. 

The same idea can be used in deriving a "Lerat method" for solving a 
general O.D.E. system like (2.18). Precisely one gets 


n+1 At 2 n+1 _ n n , , nS At ..2 n 

z + g — =— M z = z - AtMz + (1 + g) — =— M z . 


(2.31) 


The stability properties of this method are easily investigated by a normal 
mode analysis. The characteristic root of the method is 

,2 


P = 


1 - X + (1 + g) j- 


(2.32) 


1 + 8 j- 


where X = Aty and y stands for an eigenvalue of M. If 6 = ~ V 2 > a 
cancellation occurs and the characteristic root is 


P = 


ill 

“I 


(2.33) 
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hence |p | _< 1 iff >_ 0. It follows that for 3 = - 1/2 the Lerat 

discretization (2.31) of (2.18) is unconditionally stable when M is the 

matrix associated to the Chebyshev collocation spatial discretization. The 
expression (2.33) of the characteristic root shows that the Lerat scheme 

coincides with the Beam and Warming scheme on the model equation (2.1), i.e., 
whenever the coefficient matrix A is constant, so that it commutes with 
differentiation. The two methods differ for the conservation law (2.29), as 

it can be seen by comparing (2.15) and (2.30). 

If 3 * " V2 9 P has a singularity at X = ±/2/ 1 3 1 hence the scheme is 
only conditionally stable. Since the eigenvalues of the pseudospectral 
Chebyshev derivative are uniformly bounded away from the origin ([5], Sect. 
2), stability is guaranteed if At is chosen sufficiently large. 

Finally, let us observe that |p| 1 for all 3 _< “ V2 if X is 

imaginary. Hence, all the Lerat schemes are unconditionally stable when used 
with a Fourier method in space. 

Remark : The previous analysis shows that unconditional stability is 

achieved (at least for 3 = “ V2 ) is the stabilizer terra in (2.31) is built 
up by the square of the matrix of the pseudospectral derivative including the 
boundary conditions . Instead, one could think of a stabilizer term which 
simply involves the second derivative operator (in analogy to the finite 
difference case, see [3]) with no boundary condition or, say, Dirichlet 
boundary conditions. The resulting linear system would be of classical 

elliptic type, for which efficient solution techniques are available. 
Unfortunately, such a method turns out to be unstable. ||| 

Finally, we recall that the scheme (2.31) can be transformed into an 


equivalent scheme in terms of the physical unknowns using again the 
transformation (2.10). 
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3. EXTENSION TO NONLINEAR SYSTEMS 

Assume now that the system to be solved is 

+ A(w)w^ + g(w) =0, (3.1) 

where g(w) is a vector and A(w) is a 2x2 matrix which can be 

diagonalized as 

A(w) = T(w)A(w)T *(w) (3.2) 

with A(w) = diag{X,y}, X, y £ 1R, Xy < 0. The treatment of the boundary 

conditions described in the previous section can be applied here, provided the 
linear combinations of the physical equations at the boundary points involve 
as coefficients the entries of the matrix T ^(w 11 ). 

The Beam and Warming method was used to simulate a 1 hour transient of 
methane gas in a pipe of length 150 km. and diameter .5m. The equations of 
motion are (see, e.g., [1]). 


P t + (pu) x = 0 


(pu) t + 


(p + pu 


) + fpu u 

X 


= 0 


(3.3) 


where the Moody friction factor f is equal to .01. The boundary conditions 
simulate the packing of the pipe (i.e., pu is constant at the outflow and 
increasing at the inflow). 

Table I contains the values of the computed inflow pressure for different 
values of N and At. The exact value is 74.8080... The time step At has 
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the form At = aAt*, where At* is the stability limit of the Modified Euler 
method for the same problem. The CPU times (in hours) on the Honeywell 6040 
at the University of Pavia are also reported. The columns on the right 
contain the computed values and corresponding CPU times produced by the 
Modified Euler method, using At = ,9At*. 

These results show that the Beam and Warming method considered in this 
paper is unconditionally stable and more convenient than an explicit method of 
the same order for attaining accuracies of industrial interest. Further 
results will appear elsewhere. 
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TABLE I 

Comparison between Beam and Warming and Modi 
with the Chebyshev collocation method in space 


Implicit 


N 


a 


Pinf low 


CPU time 


8 


1 . 

10 . 

25. 

50. 

100 . 

200 . 

300. 


74.8084 

74.806 

74.797 

74.778 

74.509 

74.103 

73.071 


83 E - 2 
88 E - 3 
37 E - 3 
19 E - 3 
11 E - 3 
87 E - 4 
57 E - 4 


16 


1 . 

10 . 

25. 

50. 

100 . 

200 . 

300. 


74.8080 

74.806 

74.798 

74.752 

74.651 

74.274 

73.458 


32 


1 . 

10 . 

25. 

50. 

100 . 

200 . 

500. 

750. 


74.8080 

74.8080 

74.807 

74.805 

74.798 

74.752 

74.460 

74.360 


38 E - 1 

39 E - 2 

16 E - 2 
86 E - 3 
51 E - 3 
28 E - 3 
21 E - 3 

79 E0 
83 E - 1 
33 E - 1 

17 E - 1 
89 E - 2 
45 E - 2 
19 E - 2 
13 E - 2 


Euler time discretizations 


. Exact Pi n fi ow = 74.8080... 


Explicit 

Pinf low 

CPU time 

74.8076 

.58 E - 2 

74.8080 

.12 E - 1 


74.8080 


98 E -1 
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